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Two-Phase  Algorithm  for  Nonlinear  Constraint  Problems 


J.  B.  Rosen 

Computer  Science  Department 
University  of  Minnesota 


Abstract 

An  algorithm  Is  described  which  solves  the  general  nonlinear  programming 
problem  with  nonlinear  constraints.  The  computational  Implementation  Is 
based  on  the  Iterative  use  of  any  available  package  which  solves  the  linearly 
constrained  problem  with  a nonlinear  objective  function.  Large,  sparse 
problems  with  nonlinear  constraints  can  be  solved  by  this  algorithm,  provided 
a suitable  linear  constraint  package  Is  used. 

The  algorithm  consists  of  two  phases.  The  first  (Phase  I)  uses  an 
external  squared  penalty  function  to  find  a point  x^,  close  to  a local 
minimum.  Starting  with  x^,  the  algorithm  then  solves  a sequence  of  linearly 
constrained  problems  (Phase  II) . Selected  nonlinear  constraints  are  linearized 
for  each  such  Phase  II  Iteration.  With  suitable  assumptions,  convergence 
from  any  Initial  point,  with  quadratic  convergence  in  Phase  II,  is  shown. 

The  practical  Implementation  of  this  algorithm  is  described,  and  its 
potential  application  to  a model  for  the  assessment  of  energy  alternatives  Is 
discussed  briefly. 
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1.  Introduction 

Nonlinear  programming  Is  now  being  used  effectively  co  solve  relatively 
large  computer  models  for  the  assessment  of  various  energy  alternatives.  In 
a recent  study  of  energy-economy  Interactions  over  a 75  year  period  [4],  a 
model  was  developed  which  consisted  of  approximately  350  rows,  600  columns, 
and  an  objective  function  in  which  80  variables  appear  in  a nonlinear 
manner.  The  coefficient  matrix  was  somewhat  less  than  1%  dense.  A variety 
of  cases  were  run,  where  each  such  run  required  30  to  90  seconds  of  computer 
time  on  an  IBM  370/168.  The  nonlinear  aspects  of  this  model  appear  to  be 
crucial  to  a realistic  representation  of  the  economy-wide  production  function 
and  the  consumption-investment  tradeoffs.  It  also  seems  that  an  attempt  to 
represent  this  system  in  terms  of  an  entirely  linear  model  would  liave  led  to 
a larger  and  more  cumbersome  problem. 

The  solutions  were  obtained  using  the  package  MINOS  [6 >7],  which  can  be 
used  for  large,  sparse  problems  with  a nonlinear  objective  function,  but 
linear  constraints  only.  The  original  model  formulation  contained  25  non- 
linear constraints.  Fortunately,  for  the  purposes  of  the  study,  it  was 
known  that  they  were  equality  constraints,  and  each  such  nonlinear  equation 
was  used  to  eliminate  a variable  from  the  original  objective  function. 

However,  this  technique  will  not  generally  be  valid.  Therefore  an  efficient 
computer  program  capable  of  solving  large,  sparse  problems  with  nonlinear 
inequality  constraints  is  definitely  needed  for  this  class  of  problems, 
as  well  as  many  others. 

The  work  described  in  this  report  was  largely  motivated  by  this  need. 

An  Important  objective  was  to  develop  a method  for  solving  problems  with 

nonlinear  constraints,  taking  full  advantage  of  available  computer  packages 

for  solving  the  linearly  constrained  problem.  This  is  accomplished  by  the 
A 

Two-Phase  algorithm  , which  uses  such  a package  Iteratively  to  solve  a sequence 
of  linearly  constrained  problems.  During  each  iteration  die  package  is 
treated  as  a "black  box",  with  only  the  usual  Information  (linear  constraint 
matrix  coefficients,  objective  function  and  gradient  subroutines,  and  an 
initial  point)  being  supplied  to  it.  Thus  the  most  suitable  existing  package 
can  be  used,  and  if  an  improved  package  for  linearly  constrained  problems 
becomes  available,  it  can  replace  the  existing  package  with  a minimum  of 
effort. 
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* A preliminary  version  of  the  Two-Phase  algorithm  was  presented  earlier  [14]. 


1 


The  general  problem  considered  Is  the  following.  Let  ScR  denote  a 
convex  polyhedron  defined  by  the  linear  consiiralnts  and  bounds: 


Ax  ■ b 
b < X < b 


where  the  vectors  b,  b,  b and  the  elements  of  the  matrix  A are  specified. 
To  simplify  the  presentation,  It  Is  assumed  that  all  linear  Inequalities 
have  been  converted  to  equations  In  the  standard  manner.  Also,  some  (or 


all)  variables  may  be  unbounded  (b^ 


-»  or  b^  = “ , or  both).  A set  of 


q+1  functions  ()i^(x),  i*0,l,...,q,  are  defined  on  S and  (at  least  for  the 
validity  of  convergence  proofs)  are  assumed  to  have  continuous  second  deriva- 
tives for  xeS.  The  function  <^q(x)  Is  the  objective  function,  and  the 

1*1 are  the  nonlinear  constraint  functions.  The  domain  V cS  is  now 

given  by 


X e S 

4>^(x)  <0,  1*1, 


It  Is  assumed  that  V is  nonempty.  Again  to  simplify  matters,  it  Is  ssumed 
that  all  nonlinear  constraints  are  inequalities.  Nonlinear  equality  constraints 
are  handled  as  a special  case  with  no  difficulty,  as  discussed  later.  The 
general  problem  to  be  solved  is  then 

(PN)  min  (}ig(x) 

X e V 


The  theoretical  basis  for  the  solution  of  (PN)  by  the  Two-Phase  algo- 
rithm is  given  in  the  next  two  sections.  It  is  shown  there  that  with  suitable 
assumptions,  the  algorithm  will  give  global  convergence  to  at  least  a local 
minimum  x*  of  (PN).  From  an  arbitrary  initial  point  x , a point  x close 
to  X*,  is  obtained  in  Phase  I by  minimizing  an  external  squared  penalty 
function  subject  to  xeS.  An  estimate  is  also  obtained  of  the  optimal 
vector  of  multipliers  X*,  corresponding  to  x*. 

Starting  with  x^  and  X^,  the  Phase  IT  iterations  give  a sequence  of 
vectors  {x  } and  (X  }.  In  the  k*"*'  iteration  the  nonlinear  constraints  are 
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linearized  about  x and  the  minimum  is  found  to  a linearly  constrained 

problem  with  a modified  Lagrangian  as  the  objective  function.  The  solution 

1,+]^ 

to  this  problem  gives  x and  a . The  vectors  thus  obtained  satisfy 
k+1  k 

the  relation  x c M(x  ),  wViere  M is  a point-to-set  mapping.  It  is  shown 

that  a fixed  point  of  M is  a Kuhn-Tucker  point  of  (PN).  Furthermore,  with 

Ic  k 

some  additional  assumptions  it  is  shown  that  the  sequence  {x  ,X  } converges 
to  (x*,X*)  at  a quadratic  rate.  This  theoretical  justification  for  the  Two- 
Phase  algorithm  draws  heavily  on  earlier  work  on  external  penalty  functions 
[2],  and  the  iterative  linearization  of  the  nonlinear  constraints  [9,12]. 

For  a recent  discussion  of  methods  using  constraint  linearization,  and 
comparison  with  augmented  Lagrangian  methods,  see  [8j . 

The  Two-Phase  algorithm  can  be  implemented  using  any  computer  package 
which  solves  the  linearly  constrained  problem  with  a nonlinear  objective 
function.  Several  suitable  packages  for  linearly  constrained  problems  are 
now  available  [ 1,3,7,13  ].  The  Two-Phase  algorithm  has  already  been  imple- 
mented using  the  GPM  package,  and  is  being  implemented  using  MINOS. 

The  capabilities  of  the  resulting  computer  program  lor  solving  (PN) 
will  be  determined  largely  by  those  of  the  package  used.  Thus  the  size 
and  structure  of  problems  which  can  be  solved,  and  the  efficiency  with 
which  they  are  solved,  will  depend  primarily  on  the  package  used  in  the 
implementation.  The  size  and  structure  of  the  problem  is  essentially 
determined  by  that  of  the  original  linearly  constrained  domain  S,  since  in 
general,  the  additional  linearized  constraints  added  in  Phase  II  will  not 
increase  the  problem  size  significantly.  For  example,  if  the  matrix  A 
which  defines  S is  large  and  sparse,  then  the  linearly  constrained  problems 
solved  in  both  Phase  I and  Phase  II  will  also  be  large  and  sparse,  and  will 
therefore  be  solved  most  efficiently  by  a package  designed  for  such  problems. 

The  initial  implementation  of  the  Two-Phase  algorithm  has  been  completed 
using  the  GPM  package  which  is  suitable  for  relatively  small,  dense  problems 
consisting  of  no  more  than  40  variables  and  80  linear  inequality  constraints 
This  implementation  is  called  GPM/NIC.  It  has  the  advant.-ige  that  it 
is  convenient  to  use  for  relatively  small  problems,  and  It  therefore  served 
as  an  excellent  vehicle  for  developing  and  testing  the  practical  implementation 
of  the  Two-Phase  algorithm.  Using  this  package  It  was  convenient  to  experiment 
with  a number  of  possible  modifications  In  order  to  improve  the  efficiency 
of  the  resulting  computer  program. 

A variety  of  relatively  small  test  problems  have  been  run  using  GPM/NLC. 
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These  problems  all  consisted  of  15  or  fewer  variables,  and  a maximum  of  20 
constraints,  plus  bounds.  All  problems  contained  nonlinear  Inequality 
constraints,  both  convex  and  nonconvex.  Some  also  had  nonlinear  equality 
and  linear  constraints.  A number  of  these  problems  have  been  used  previously 
as  test  problems  for  other  methods  [16] . In  all  cases  the  performance  by  GPM/NLC 
was  comparable,  and  often  better,  than  that  of  other  methods  on  the  same 
problems.  This  performance  depends  on  the  proper  choice  of  the  penalty 
parameter  p In  Phase  I,  and  on  the  quadratic  convergence  rate  achieved  In 
Phase  II.  Normally  convergence  is  attained  with  no  more  than  five  iterations 
in  Phase  II,  with  the  final  two  of  these  relatively  fast,  since  the  corres- 
ponding change  In  x Is  small. 

A preliminary  implementation  of  the  Two-Phase  algorithm  using  MINOS  has 
also  been  tested.  This  implementation  Is  called  MINOS/NLC.  These  tests  have 
used  a modified  version  of  the  energy  alternative  model  discussed  above  witli 
350  linear  constraints.  In  addition  a total  of  30  nonlinear  inequality  con- 
straints have  been  added,  only  some  of  which  are  active  at  the  optimum. 

Based  on  these  preliminary  tests  it  is  estimated  that  nonlinear  constraint 
problems  of  this  size  can  be  solved  in  3 to  5 minutes  on  an  IBM  370/168,  or 
comparable  machine.  Details  of  the  computational  results  obtained  with 
GPM/NLC  and  MINOS/NLC  will  be  given  elsewhere  [15]. 
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2.  Linearization  of  Nonlinear  Constraints 


We  first  consider  the  linearly  constrained  minimizatLon  problem 
(PS)  min  •J'qCx) 


Since  S is  bounded,  there  is  a point  x*  e S at  which  '('qCx)  attains  its  global 
minimum;  that  is,  ijj^Cx*)  ^ 4>q(x)  for  all  xeS.  It  follows  directly  from  the 
previous  assumptions  that  a global  minimum  of  (PN)  is  attained  at  some  point 
X*  e V,  which  may  not,  however,  be  unique.  Since  VcS,  we  have  (f^Cx*)  _<  cJi^Cx*) 
A special  situation  occurs  if  x*  e V,  in  which  case  we  have  x*  = x*  (all  the 
nonlinear  constraints  are  redundant). 

While  we  would  like  to  determine  a global  minimum  of  (PN) , we  will 
frequently  only  find  a local  minimum.  Provided  that  the  4*^  satisfy  some 
regularity  condition,  a local  minimum  will  always  be  a Kuhn-Tucker  point 
(i.e.,  the  first-order  Kuhn-Tucker  optimality  conditions  will  be  satisfied). 

Now  given  any  point  y e S,  we  consider  the  1 inearized  approximations  to 
the  constraints  <{i^(x)  ^ 0.  That  is,  we  replace  them  with 

h^(y;x)  H (<i^(y)  + (x-y)^74)^(y)  £ 0,  i = l,...,q 
This  gives  us  a polyhedral  set  W(y)c:  S defined  by 


p I h^(y;x)  _<  0,  i=l qj 

Note  that  if,  say  4’^  (x)  = 0,  we  use  the  corresponding  linear  equality 
hj(y;x)  = 0 in  W(y).  For  an  arbitrary  y s S,  W(y)  could  be  empty. 

However,  W(y)  is  not  empty  for  any  y e V,  since  yeV  implies  yt:W(y).  It  is 
convenient  to  consider  W as  a point-to-set  mapping  W:S  -*■  S.  We  will  make  an 
additional  assumption  later  about  the  functions  4^^,  to  insure  the  continuity 
of  the  mapping  W. 

The  linearized  constraints  give  us  a problem  which  is  closely  related 
to  (PN).  Given  any  y e S we  consider  the  problem 

(PL)  min  4'(y;x) 
xeW(y) 
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W-T  ' ' 

where  'l'(y;x)  = ^ o(y;x),  and  a(y;x)  satisfies  the  relations  o(y;y)  ■ 0 

and  7o(y;y)  • 0.  Note  that  this  problem  consists  of  linear  constraints  and 
bounds  only,  with  the  nonlinearity  confined  to  the  objective  function. 

As  with  (PN)  we  would  like  to  find  a global  minimum  of  (PL),  but  we  will 
often  find  only  a local  minimum.  We  therefore  Interpret  an  optimal  solution 
to  (PL)  as  any  point  xeW(y)  at  which  'i'(y;x)  has  a local  (constrained  or 
unconstrained)  minimum.  In  general,  the  word  minimum  will  mean  a local 
minimum,  unless  global  minimum  Is  explicitly  stated. 

The  problem  (PL)  will  always  have  an  optimal  solution,  provided  that 
W(y)  Is  not  empty.  The  problem  (PL)  therefore  defines  another  point-to-set 
mapping  M:S  -*■  S.  Specifically,  M(y)‘=W(y)  is  the  point  (or  set  of  points)  at 
which  'f(y;x)  attains  Its  minimum  for  xeW(y).  That  is 

M(y)  5 arg  min  4'(y;x) 
xcW(y) 

Our  main  Interest  Is  In  a fixed  point  of  this  mapping,  that  Is  a point 
y such  that  yeM(y). 

Theorem  1 

A fixed  point  of  the  mapping  M is  a Kuhn-Tucker  point  of  (PN) . 

Proof:  Let  y be  a fixed  point  of  M.  Then  we  have  yeW(y).  We  first  show 
that  y e V.  If  y ^ V,  then  for  at  least  one  i,  say  i = k,  we  must  have  <l'j^(y)  ^ 0- 
But  then  we  would  have  hj^(y;y)  = so  that  y^W(y). 

Since  y is  an  optimal  point  for  the  problem  (PL)  with  all  constraints 
linear,  it  follows  that  the  first-order  Kuhn-Tucker  optimality  conditions 
are  satisfied  at  y.  These  conditions  Involve  the  functions  hj^(y;x)  and  their 
gradients  Vh^(y;x)  evaluated  at  x=y,  for  all  active  linearized  constraints 
at  y,  that  is,  those  for  which  h^(y,y)=0.  The  active  linear  constraints  and 
bounds,  as  well  as  the  objective  function  gradient  Vf(y;y),  also  enter  into 
these  conditions.  But  these  are  just  the  Kuhn-Tucker  optimality  conditions 
for  the  problem  (PN)  at  the  point  X“y,  since  we  have  h^(y;y)  “ 4'^(y)» 

^h^(y;y)  = 7<l(^(y),  and  V4'(y;y)  « V(jiQ(y).B 

It  should  be  noted  that  we  also  have  the  same  value  for  the  objective 
functions  of  (PL)  and  (PN)  at  a fixed  point  y,  since  '4'(y;y)  = liiQ(y). 

We  now  consider  using  the  mapping  M as  the  essential  part  of  an  algorithm 
to  solve  (PN).  For  this  purpose  we  desire  an  algoritlim  which  generates  a 

sequence  of  points  converging  to  a fixed  point  of  M.  We  consider  the  following. 
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Algorithm  A 

0 k ' 

Starting  with  an  initial  point  x eS,  obrain  the  sequence  of  points  {x  },  i 

generated  by  (PL).  For  this  sequence  we  have 

k+1  k < 

(1)  X ^ eMCx*^)  , k = 0,1,...  j 

i 

The  convergence  of  this  algorithm  to  a fixed  point  will,  in  general,  1 

depend  on  the  use  of  an  appropriate  term  a(y;x)  in  the  objective  function. 

It  is,  however,  of  some  interest  to  show  that  for  a suitably  restricted  class 
of  problems  this  alg'oritiim  converges  to  a fixed  point  of  M,  starting  with 
any  x^  e V,  even  for  a = 0. 

Assume  that  the  are  concave  on  S,  and  that  they  satisfy  a regularity 
condition  which  ensures  that  the  mapping  W is  continuous  [5,10].  Note  that 
for  concave  the  feasible  set  V will  not,  in  general,  be  convex,  and  in 
fact  can  be  called  reverse-convex  [5].  For  an  earlier  discussion  of  this 
case  see  [11]. 

0 k 

Lemma  1 Starting  with  any  x c V,  generate  a sequence  {x  } by  means  of 

k k+1 

Algorithm  A,  with  a = 0.  Then  for  k = 0,J,...,  we  have  x eV  and  (x  ) ^ 
k 

(|!q(x  ).  Furthermore,  there  will  be  a convergent  subsequence,  which  converges 
to  a fixed  point  of  M. 

Proof:  For  concave  we  have  if.(x)  jc  h^(v,x)  for  any  x.yeS.  Then  for  any 
X e W(y)  we  have  hj(y;x)  <0,  i=  1, . , . ,q,  which  implies  -t.  (x)  £ 0,  1=1,...  ,q. 

Thus  yeW(y)c:v,  for  any  y c V. 

k. 

Since  V<=S  is  compact,  the  sequence  {x  } has  an  accumulation  point,  say 

k 

X*  c V.  Then  because  of  monotontcity  and  continuity,  ((-^(x*)  < ^qIx  ) for  all 

k.  Now  suppose  x*  is  not  a fixed  point  of  M.  Then  there  must  exist  a point 

xcW(x*),  with  X close  to  x*,  such  that  (l)„(x)  < 3y  the  continuity 

J k ^ 

of  the  mapping  M there  then  would  be  a point  x of  the  sequence,  sufficiently 
close  to  X*  such  that  (fi^Cx^^^)  < (fi^Cx*),  a contradiction. Bf 

The  well-known  difficulty  with  Algorithm  A for  a = 0,  is  that  it  will 
not,  in  general,  solve  (PN)  for  the  standard  convex  problem  (ail  convex), 
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where  any  Kuhn-Tucker  point  Is  a global  minimum.  This  can  be  easily  seen 

2 

by  applying  it  to  the  problem  in  R : 


min  . 


■*l"*2 


0 £ Xj^,  X2  <.  2 

2 2 
+ x^  < 2 


with  the  optimum  solution  x*  ■ x*  = 1.  For  any  starting  point  not  on  the 
line  Xj^  •=  X2»  the  sequence  of  points  generated  will  all  lie  on  the  boundaries 
of  the  square. 

This  difficulty  can  be  overcome  by  a suitable  choice  for  o.  Replacing 
<|)q  bv  the  Lagranglan  [12],  or  a modified  Lagrangian  [ 9 ] will  ensure  con- 
vergence (in  fact,  at  a quadratic  rate)  provided  the  starting  point  is  close 
enough  to  a local  minimum.  The  modified  Lagranglan  will  be  used  here.  This 
is  obtained  by  choosing 

(2)  o(y;x)  - f ^^(y)  [i[>^(x)  - h (y;x)] 
i®! 

where  the  Lagrange  multipliers  X are  determined  in  an  appropriate  manner. 

^ k k k 

More  specifically,  the  solution  of  (PL)  with  y“x  , and  X.(x  ) = X , gives  a 

k+1  k k+1  ^ ^ 

point  X e M(x  ) and  also  a set  of  multipliers  X^  i l = l,...,q,  cor- 
responding to  the  linearized  constraints.  Only  those  multipliers  which 
correspond  to  an  active  linearized  constraint  can  be  positive,  and  therefore 
contribute  to  a. 

We  consider  a local  minimum  x*  of  (PN)  and  let  X*  ^ 0 be  the  corresponding 
vector  of  multipliers.  We  make  a regularity  assumption  about  the  constraints 
at  X*  so  that  x*  is  a Kuhn-Tucker  point,  and  also  assume  that  second-order 
sufficiency  conditions  hold  there.  Some  additional  assumptions  are  needed 
to  insure  that  certain  matrices  (Involving  the  Lagrangian  and  its  derivatives) 
are  bounded  in  the  neighborhood  of  (x*,X*). 

k th 

Lemma  2.  Apply  Algorithm  A,  with  a as  given  by  (2)  and  y=x  at  the  k*^ 

iteration,  to  generate  a sequence  of  points  {x  } and  corresponding  vector 

k 11 

of  multipliers  {X  },  k=2,3 starting  with  specified  vectors  x and  X . 

Then  if 


the  sequence  c will  converge  to  (x*,X*)  at  a quadratic  rate.  The 

constants  o and  a depend  only  on  the  norm  of  certain  matrices  evaluated  in 
the  neighborhood  of  (x*,X>'). 

More  specifically,  we  have 

-.k 


(4) 

(5) 

(6) 


< 2(|)^  ||f(:5^-) 


I Az 

I k 


k+1  i 


< a Az 


ki  ,2 


2 12 

< -ib 

— OL  2 


for  k=l,2,...  . In  these  relations  wa  simplify  notation  by  using  z (x.X) 
ri'^Q  k k*^l  k 

£ R Az  r z - z . The  vector  function  f(z)  Is  given  by 

£(z)  SR"-"’ 

\w  (x,A^ 

where  A(x,X)  = IX^<|)^(x)  is  the  Lagrangian  for  (?N)  , A'  is  its  gradient 

with  respect  to  x,  and  w^(x,X)  = (X  (j)  (x)  , . . . , X (})  (x)), 

i 1 q q 

The  detailed  assumptions  and  proof  have  beer,  given  by  Robinson  [ 9 ) , 
and  will  not  be  repeated  here.  The  results  hold  with  ] j ' | | representing 
any  consistent  norm.  To  simplify  the  presentation  here,  it  has  been  assuried 
that  only  nonlinear  constraints  are  active  at  x*.  If  any  of  the  original 
linear  constraints  or  bounds  are  active,  the  gradient  V((i^  must  be  replaced 
by  its  projection  in  the  appropriate  subspace. 

It  should  be  noted  that  provided  X^  > 0 whenever  !j>^(x)  > 0,  the  point 
(y*,X*)  is  a Kuhn-Tucker  point  if,  and  only  if,  f(x*,X*)  = 0. 

Thus  from  a good  starting  point  Algorithm  A gives  rapid  convergence  to 
the  closest  local  minimum.  However,  this  still  does  not  give  us  a satisfactory 
method  for  getting  a Kuhn-Tucker  point  from  a more  or  less  arbitrary  given 
point  X . To  achieve  this,  we  must  first  generate  a good  starting  point 
X and  a corresponding  multiplier  X'^.  For  this  purpose  we  need  a method 
which  generates  a point  close  to  a local  minimum  of  (PN),  from  an  initial 
point  x*^  which  may  lie  anywhere  in  S (the  original  linearly  constrained 
domain).  Since  the  linearization  of  the  about  such  an  x*^  may  give  a 
very  poor  approximation  to  the  nonlinear  constraints,  we  do  not  want  to  use 
(PL)  initially. 
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3 . External  Squared  Penalty 


The  most  satisfactory  answer  appears  to  be  the  use  of  an  external 
squared  penalty  for  the  nonlinear  constraints.  That  Is,  we  solve  the  linearly 
constrained  problem 


(PI) 


min  H'(ij;x) 
xeS 


where 


5 ? 

1-1 


(x)r 


The  for  > 0,  and  are 


j 


zero  otherwise.  Nonlinear  equality  constraints  are  Included  by  using  4’ 

itself  (rather  than  4i^)  if  4> . (x)  = 0. 

^ ^ 1 
The  penalty  y must  be  chosen  large  enough  to  Insure  that  x is  sufficiently 

close  to  X*,  but  no  larger  than  necessary,  since  a large  value  of  y can 

greatly  Increase  the  computation  time  of  (PI).  The  solution  of  (PI)  also 

gives  an  estimate  of  the  optimal  vector  of  multipliers  A*. 

The  use  of  the  external  squared  penalty  to  solve  (PN)  has  been  thoroughly 

investigated  [ 2 ] . By  choosing  a sequence  of  values  {y  }-><»,  it  is  shown 

that  the  corresponding  sequence  of  points  {x  },  given  by  (PI)  converges  to 

a local  minimum  of  (PN) . The  conditions  on  (PN)  for  convergence  from  an 

arbitrary  initial  point  appear  to  be  the  least  restrictive  of  any  globally 

convergent  algorithm.  For  example,  a condition  which  will  ensure  global 

convergence  to  at  least  a local  minimum  of  (PN)  is  that  the  function 


i=l 


+ 2 

[i{i^(x)]  be  quaslconvex  on  S.  We  assume  some  such  condition  to  hold. 

However,  instead  of  solving  a sequence  of  problems  (PI)  with  increasing 
values  of  y,  we  pick  a single  (relatively  small)  value  of  y and  solve  the 
corresponding  problem  (PI)  to  get  a point  x^  = x(y) . We  call  this  Phase  I. 

In  general,  x^  will  violate  some  of  the  nonlinear  constraints.  We 
estimate  A*  by  computing  the  vector  of  multipliers  A^  = A(y),  given  by 

(7)  I^(u)  = y<(i^(x(y))  >0,  i-l,...,q. 

Once  we  have  an  x^  and  A^  obtained  in  this  way  we  can  use  the  previous 
algorithm  to  obtain  rapid  convergence  to  a local  minimum.  We  call  this  Phase 
II.  The  combination  of  Phase  I and  Phase  II  leads  directly  to: 
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Algorithm  B 

1.  Start  with  x^eS. 

2.  Phase  I:  Solve  (PI)  with  a specified  u,  to  get  x and  A as 

given  by  (7).  Set  x^,  a'^  A^  and  k l. 

k k 

3.  Phase  II:  Given  x and  A , solve  (PL)  with  a as  given  by  (2) 

k k+1  k+1 

and  y = x . This  gives  x and  A 21  0. 

/ rr  I I / \k+l>  / k .k.  II  k+1  ,k+l 

4.  If  I I (x  ,A  ) - (x  ,A  )||  ^e,  output  x and  A . Stop. 

CP  k k+1  , k k+1  , , . 1 j ^ 1 

5.  Set  x-»-x  jA-^A  ,k-*-  k+1,  and  return  to  3. 


We  will  now  prove  convergence  of  this  algorithm  to  a local  minimum  of  x*  of 
(PN)  from  any  point  x^ c S.  This  proof  is  based  on  the  assumption  that  starting 
from  any  x e S the  external  squared  penalty  method  (PI)  converges  to  a local 


minimum  x*,  as  p 
be  denoted  by  A*. 
2 are  valid. 


«>.  The  corresponding  optimal  vector  of  multipliers  will 
We  also  require  that  the  assumptions  made  prior  to  Lemma 


Theorem  2 

There  is  a y > 0,  such  that  if  we  use  any  p p in  Phase  I,  and  starting 

0 k k 

with  X eS  generate  a sequence  of  vectors  {x  ,A  },  k=l,2,...,  by  Algorithm 

B,  this  sequence  will  converge  to  {x*,A*}  with  a quadratic  convergence  rate. 

1 k+1  k+1 

For  any  e Phase  II  will  terminate  with  vectors  x and  A which 

satisfy 

(8) 


k+1  ,k+l 


(x  , A ) - (x*,  A*)  I I £ e 


Proof;  For  every  u > 0,  (PI)  determines  an  optimal  vector  x(p)  and  the 
corresponding  A(p)  given  by  (7).  By  assumption  lim  x(u)  = x*.  As  shown  in 

[2],  we  also  have  that  lim  A(p)  = A*.  Thus  we  can  choose  p so  that  with 


x^  = x(p)  and  A^  = A(p),  we  satisfy  the  condition  O'*  of  Lemma  2.  By  Lemma 
k k 

2,  the  sequence  {x  ,A  } converges  quadratically  to  (x*,A*).  At  termination 
we  have  | ]Az  1 | ^ g ^ • Then  using  (5),  it  follows  that 


, k+1  k+1.  . . . .. 

(x  ,A  ) - (x*,A*) 


j-k+1 


I 1 1 ^ 1 1 E 

j=*0 
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4 . Practical  Implementation 


r 

it 

\ 


I 
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In  the  preceding  sections  the  theoretical  basis  for  the  Two-Phase 
algorithm  has  been  presented.  We  now  consider  the  practical  Implementation 
of  this  algorithm  to  give  an  efficient  computer  program  for  finding  at 
least  a local  minimum  of  the  problem  (PN) . In  addition  to  efficiency,  It 
Is  desired  to  accomplish  two  related  objectives  with  this  Implementation. 

First,  the  Implementation  should  be  able  to  use  any  existing  package 
capable  of  solving  the  linearly  constrained  nonlinear  programming  problem. 
This  allows  the  user  to  take  full  advantage  of  the  structure  of  the  linear 
constraint  matrix  A (relatively  small  and  dense,  large  and  sparse,  etc.) 
by  using  a suitable  package.  It  also  permits  the  relatively  easy  substitu- 
tion of  Improved  packages  as  they  become  available.  Some  existing  packages 
(MINOS  in  particular)  require  very  little  from  the  user  in  the  way  of 
tolerance  or  parameter  selection.  This  benefit  will  then  carry  over  directly 
to  the  nonlinear  constraint  problem,  when  such  a package  is  used. 

A second  important  objective  is  a robust  program,  in  the  sense  that 
it  will  find  at  least  a local  minimum  for  a variety  of  problems  typical  of 
those  found  in  practical  applications.  Typically  the  objective  function  and 
at  least  some  of  the  constraints  will  not  be  convex,  only  some  of  the  con- 
straints will  be  active  at  the  minimum,  and  the  initial  point  may  not  be 
feasible.  Another  practical  difficulty  will  arise  if  the  set  S is  feasible, 
but  there  is  no  feasible  solution  to  (PN) , that  is,  V is  empty.  This  could 
happen  as  a result  of  conditions  which  are  too  stringent,  unrealistic 
assumptions,  or  an  error  in  formulation.  In  any  event,  a practical  imple- 
mentation cannot  assume  that  a feasible  solution  to  V exists,  but  must  have 
provision  to  recognize  this  difficulty,  if  it  occurs,  with  a minimum  of 
wasted  comt>uter  time. 

The  modifications  to  Algorithm  B which  are  described  in  this  section  are 
intended  to  achieve  these  practical  objectives.  In  describing  the  modified 
algorithm  (to  be  called  Algorithm  C)  we  will  assume  that  the  package  used 
to  solve  the  linearly  constrained  problems  has  certain  desired  properties. 
These  are: 

1.  The  ability  to  find  a feasible  point  with  respect  to  the  specified 
linear  constraints  and  bounds,  if  such  a point  exists.  If  no  such  point 


i 


fcxlHts,  a "no  feaMltio  solatlop"  exit  from  the  package  Ib  activated. 

2.  The  package  will  iicrmaliy  find  at  least  a locf.l  .-ninimun:  for  any 
differentiable  objective  function,  subject  to  linear  constraints  and 

bounds.  The  optimality  test  rei’utres  that  the  first  order  Kuhn-Tucker  ' 

conditions  be  satisfied  within  <'i  certain  tolerance.  This  tolerance  will  i 

depend  on  certain  computed  quantities  and  may  also  depenrt  on  a user  specified  ' 

tolerance.  When  the  optimality  test  is  satisfied  the  optimal  vector  x, 

.ai'.d  corresponding  vector  X of  multipliers,  is  available.  If  the  optimality  j 

test  is  not  satisfied  after  a specified  maxiraam  number  of  iterations,  the 
calculation  Is  stopped,  and  the  current  vectors  x and  X are  available. 

3.  The  package  must  be  able  to  accept  an  initial  vector  x^,  and  if 
Is  a local  minimam  It  should  recognize  this  fact  without  unnecessary 

computation.  If  is  not  a local  minimum  It  should  find  a "close"  local 
mlnimumi  that  is,  if  is  feasible  it  should  find  a sequence  of  points 
with  decreasing  function  values,  and  If  x^  is  not  feasible  it  should  first 
find  a "close"  feasible  point  and  then  continue  downhill. 

4.  The  package  requires  only  the  evaluation  of  the  objective  function 
and  its  gradient.  Subroutines  for  these  evaluations  are  normally  supplied 

by  the  user,  although  the  gradient  components  may  be  obtained  by  differencing 
In  some  cases. 

At  least  two  such  packages  are  now  available  [7,13].  We  will  refer  to 
this  package  for  linearly  censtrained  problems  as  the  IX  package. 

We  now  present  the  modified  algorithm  designed  to  achieve  the  practical 
objectives  discussed  above.  The  algorithm  will  be  given  first,  followed  by 
definitions  of  terms  not  previously  defined.  A more  detailed  discussion  of 
the  motivation  for  the  modifications  and  remarks  on  the  algorithm  are  given 
in  Section  5. 

Algorithm  C 

Input  Data:  matrix  A;  vectors  b,  b,  b;  functions  <!>  , l = 0,l,...,q, 

0 ~ ^ 

and  their  gradients;  initial  point  x ; tolerances  e,  6;  parameters  y,  k. 

Also,  any  user  specified  tolerances  and  parameters  required  by  the  ^ package. 

Q 

1.  Phase  1:  Generate  the  linearly  constrained  set  (x  ),  for  input 

1 1 

^ 


to  (PIM). 


2. 

Call  IX  package  and 

and  X as  defined  by  (9).  If 

vectors. 

and  Stop. 

3. 

o _ k 1 >k  .1 

Set  X X , X ■«-  X 

4. 

Phase  II;  Given  x 

Wjj(x^, 

for  input  to  (PLM) . 

5. 

Call  LC  package  and  i 

k+1 

,k+l 

X and 

X > 0. 

6. 

If  1 |x^^^  - x^l 1 ^ e 

Stop. 

7. 

I£  i 11 

8. 

If  k * k,  Stop. 

9. 

^ k k+1  ,k 

Set  X X , \ 

k+1  ,k+l. 


k. 


Definitions  of  terms  used  in  Algorithm  C: 

The  set  of  constraints  linearized  in  Phase  I depends  on  x^,  and  is 
given  by: 

5°  = {i  I |<(.^(x°)i  < 10^6} 

The  linearly  constrained  set  for  Phase  I is  then  defined  as  follows: 


Wj(x  ) = -jx 


X G S __  I 

h^(x®;x)  _<  0,  1 e J ' 


The  problem  solved  in  Phase  I is  given  by: 


(PIM)  min  i(n;x) 

xeWj^  (x^) 


where 


and 


y(li;x)  H (>q(x)  + 


1 ? I*: 
^ 1-1  ^ 


(x)]' 


u = 2y/10  6 
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The  solution  to  (PIM)  gives  x"^  and  a set  of  multipliers  ^ 0,  corresponding 

to  those  constraints  h.  (x*^;x)  £ 0,  IcJ^,  which  are  active  at  x^.  We  define 

^ 1 

the  corresponding  vector  X , with  components  given  by 


(9) 


X^  = max{X^,ij^^(x^)}  , l-l,...,q 


The  corresponding  sets  for  Phase  II  are  defined  as  follows: 

i (1  I -10^6  < (ji^(x^)} 


/ -= 

j-1 


Wii(x  ) = -{x 


X £ S 

h.  (x*^,x)  <0,  1 e 
1 ~ c 


The  problem  solved  In  the  kth  major  Iteration  of  Phase  II  Is  then  given  by: 

(PLM)  min 


xeWjj(x  ) 


where 


'i'’^(x)  = 4*0 (x)  + 1 - h^(x'^;x)] 


1=1 

k+1  k+1 

The  solution  to  (PLM)  gives  x and  X 2.  package  exit  from 

(PUl)  occurs  as  a result  of  the  specified  Iteration  limit,  the  vector  X 

k+1 

available  may  have  some  negative  elements.  In  that  case  we  let  X^  = 
max{0,X^}. 


15 


5.  Discussion 


Ue  conclude  this  paper  with  a more  detailed  discussion  of  the  motivation 
for  the  modifications  required  In  practical  Implementation  of  the  algorithm, 
and  also  some  additional  remarks  on  the  behavior  of  Algorithm  C based  on 
computational  experience  with  a variety  of  test  problems. 

In  order  to  achieve  computational  efficiency.  It  is  important  to  use  a 
proper  value  of  p In  Phase  I.  A large  value  of  p will  cause  the  Hessian 
matrix  of  ^ to  be  Ill-conditioned,  resulting  In  a large  increase  In  the  com- 
putational effort  (measured  in  either  function  calls  or  computer  time)  In 
Phase  I.  If  too  small  a value  Is  used,  the  point  x^  (and  corresponding  X^) 
obtained  in  Phase  I will  not  be  close  enough  to  a local  minimum  (x*,X*) 
to  ensure  convergence.  Fortunately,  there  is  a range  within  which  any  value 
of  p will  be  satisfactory.  This  is  Illustrated  in  Figure  1,  which  shows  the 
(averaged  and  smoothed)  behavior  for  some  typical  15  variable  problems  using 
the  GPM/NLC  implementation.  The  three  curves  show  the  number  of  function  calls 
required  in  Phase  I,  in  Phase  II,  and  the  total  in  both.  F.ach  function  call 
requires  the  evaluation  of  the  function  and  its  gradient  corresponding  to 

and  each  of  the  active  constraint  functions  (t>.  • The  Phase  II  curve  includes 

^ 11 

all  the  iterations  required  in  Phase  II.  As  p is  increased,  the  x and  A 
given  by  Phase  I become  closer  to  x*  and  X*,  so  that  fewer  Phase  II  iterations 
are  required.  The  curve  of  total  function  calls  is  seen  to  have  a relatively 
broad  minimum,  so  that  an  order  of  magnitude  change  in  p from  its  best  value 
would  cause  only  a small  increase  in  the  computational  effort  required.  For 
a different  type  and  size  of  problem  the  shape  of  these  curves  would  be 
similar,  but  the  coordinate  values  would,  in  general,  be  different. 

Roughly  speaking,  the  best  choice  for  p is  the  smallest  value  which  will 
produce  an  (x^,X^)  that  gives  convergence  in  Phase  II.  It  is  well  known  that 
one  often  gets  convergence  of  Newton's  method,  as  applied  to  solving  a system 
of  nonlinear  equations,  with  a starting  point  which  is  reasonably  good,  but 
does  not  satisfy  a sufficient  condition  for  convergence.  A similar  situation 
may  occur  here. 

We  take  advantage  of  this  situation  by  choosing  p just  large  enough  so 
that  at  the  end  of  say  k iterations  in  Phase  II  we  will  satisfy  all  nonlinear 
constraints  within  a specified  tolerance  6.  A satisfactory  value  of  k would 
normally  be  k = 5. 


J 
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From  (7)  we  have 


I |^(w)  I I •=  u|  l4''^(x(u))  I I 

where  , iji^)  e . Since  X(u)  is  an  approximation  to  X*,  we  have 

I I ()>''’ (x(p))  j j _<  2y/p 

where  y is  an  upper  bound  on  | |X*| |.  Now  suppose  we  choose  p = p * 2y/10^6. 
Then  at  the  end  of  Phase  I the  corresponding  x^  = x(p)  will  be  such  that 
lU^Cx  )||  £ 10^6.  From  (4)  we  get  the  approximate  relation  for  Phase  II 
iterations: 


< 2(i)^  lir(xh 


Thus  for  any  k > k = 5,  we  have 


l>'^(x*^‘^^)  I I < 2xlo^(i)^  6 < 0.0056 


Based  on  computational  experience,  the  choice  p = p = 2y/10  6 in 
Phase  I,  achieves  the  desired  convergence  in  Phase  II.  However,  it  will  not 
always  give  an  (x^,X^)  in  the  quadratic  convergence  neighborhood  of  (x*,X*). 
Therefore,  the  first  few  iterations  of  Phase  II  may  not  satisfy  (3).  This  will 
will  not  cause  any  practical  difficulty,  provided  only  that  at  least  the 
last  iteration  of  Phase  II  is  taken  in  the  quadratic  convergence  neighborhood. 

In  that  case.  Theorem  2 applies  and  the  error  bound  (8)  is  valid. 

The  nonlinear  constraint  tolerance  6 must  be  chosen  by  the  user,  since 
it  requires  knowledge  of  the  significance  of  a small  change  in  the  value  of 
a constraint  function  Instead  of  the  requirement  <Ji^(x)  ^ 0,  a constraint 

is  assumed  to  be  satisfied  if  i)>^(x)  _<  6.  In  effect,  the  user  specifies  a 
tolerance  such  that  a change  in  the  constraint  value  (or  the  equivalent,  its 
right-hand-side)  of  magnitude  6,  or  less,  is  considered  to  be  Insignificant 
from  a practical  point  of  view.  Constraints  are  scaled  where  necessary,  so 
that  the  same  tolerance  6 holds  for  each  of  them. 

This  tolerance  also  defines  the  active  constraints.  The  constraint 
is  said  to  be  active  at  a point  x,  if  1 (x)  j ^6.  In  order  to  reduce 
unnecessary  computation,  only  those  constraints  are  linearized  that  might  be 
part  of  the  active  set  at  the  optimum.  This  is  accomplished  by  only  linearizing 
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chose  constraints  which  are  either  violated,  or  come  close  to  being  active 
at  the  end  of  Phaae  I,  or  any  Phase  II  iteration.  The  Justification  for  this 
Is  Chat  If  we  knew  In  advance  which  nonlinear  constraints  were  Inactive  at 


X*,  that  Is  those  for  which  ^^(x*)  < -6,  we  could  eliminate  them  from  con- 
sideration entirely  without  changing  the  problem  (at  least  locally). 

Another  difficulty  which  can  arise  from  linearizing  a constraint 

about  X,  when  <^.(x)  Is  negative  and  relatively  large  In  magnitude,  is  shown 
1 2 2 
by  the  following  simple  two  variable  problem.  Let  <ti  = - x^  + 1 5,  0,  be 

Che  nonlinear  constraint  (outside  of  circle  Is  feasible).  Also  assume  the 

bound  Xj^  ^ 1.5.  .\ny  point  (Xj^,X2)  outside  the  circle  of  unit  radius  and 

with  Xj^  1.5,  Is  of  course  feasible.  However,  It  Is  easy  to  see  Chat 

linearizing  ♦ about  (x^,x^)  = (3,0)  gives  x^  5/3,  so  that  there  is  no 

feasible  solution  to  Che  linearized  problem. 

On  the  other  hand.  It  is  useful  to  linearize  some  constraints  even  in 

Phase  I,  If  we  wish  to  recognize  a local  minimum  (or  a point  close  to  a 

local  minimum)  when  such  a point  is  used  as  the  initial  point  x*^  in  Phase  I. 

This  difficulty  is  resolved  by  linearizing  all  those  constraints  which  are 

close  to  x*^,  more  specifically  those  for  which  ieJ*^.  Even  if  a constraint, 

^ 1 

say  iji.,  has  been  linearized  in  Phase  1 and  in  addition  is  active  at  x (with 
_ J 1 

Xj  >0),  it  is  likely  that  we  will  have  4’j  ^ 0*  view  of  this  possibility, 

the  equation  (9)  gives  the  best  estimate  for  X*  at  the  end  of  Phase  I. 

We  now  conclude  with  some  remarks  which  may  be  helpful  in  understanding 
Algorithm  C. 

1.  If  x^  is  a local  minimum,  then  all  nonlinear  constraints  active  at 

this  minimum  (i.e.,  those  with  (x^)=0)  will  be  linearized  and  included 
0 0 ^ 

in  W (x  ).  Therefore,  x will  be  an  optimal  solution  to  (PIM) , so  that 
^ 10 

it  will  give  x =x  . Thus,  if  the  initial  point  is  a local  minimum,  it 
will  be  recognized  as  such  at  the  start  of  Phase  1. 

2.  The  Initial  point  x^  need  not  be  feasible.  If  the  original  linear 
constraint  set  S has  no  feasible  solution,  this  will  be  recognized  by  the 
LC  package  at  the  start  of  Phase  I (i.e.,  a no  feasible  solution  exit 
from  LC) . 


3.  If  S Is  not  empty,  the  sequence  of  points  {x  } generated  in  Phase 
I and  II  will  all  be  in  S.  However,  even  if  x eV,  we  will  usually  not 


f 


r 

f 


Ic  *•  ^ Ic 

have  X e V,  k - 1,2, . . . ,1c.  Thus,  in  general,  | | (^  (x  ) | | > 0,  although  this 

constraint  violation  norm  will  go  to  zero  quadratically  as  given  by  (10). 

4.  The  LC  package  will  legitimately  stop  for  one  of  three  reasons: 

(a)  no  feasible  solution  exists  to  the  linear  constraints;  (b)  an  optimal 
solution  Is  obtained  within  the  LC  tolerances;  (c)  the  LC  Iteration  limit 
is  reached. 

If  (a)  occurs  it  indicates  no  feasible  solution  exists  to  (PN) , or 

other  serious  difficulty,  and  the  algorithm  terminates.  If  (b)  or  (c) 

k+1  k+1 

occur,  the  vectors  x and  X are  available  and  the  algorithm  continues. 


Successful  termination  of  the  algorithm  occurs  only  In  Step  2 or 


Step  6. 


6.  If  the  algorithm  terminates  in  Step  7 or  8,  the  vectors  x and 
k+1 

X are  available,  but  in  general  are  not  optimal. 


7.  A constraint  function  <t>.  and  its  gradient  is  only  called  in  (PIM) 

+ ^ k 

when  > 0.  Similarly  in  (PLM) , only  when  X^  > 0.  This  can  significantly 

reduce  the  amount  of  computation  required. 

8.  The  constants  10^  used  in  the  definition  of  and  J'^ , and  10^  used 
in  computing  p were  determined  as  best  so  as  to  minimize  the  computational 
effort  in  a variety  of  test  problems.  A change  in  these  values  by  a 
factor  of  10  seems  to  have  little  effect. 
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77-22  "Two-Phase  Algorithm  for  Nonlinear  Constraint  Problems" 

An  algorithm  Is  described  which  solves  the  general  nonlinear  programming 
problem  with  nonlinear  constraints.  The  computational  Implementation  Is 
based  on  the  Iterative  use  of  any  available  package  which  solves  the  linearly 
constrained  problem  with  a nonlinear  objective  function.  Large,  sparse 
problems  with  nonlinear  constraints  can  be  solved  by  this  algorithm,  provided 
a suitable  linear  constraint  package  Is  used. 

The  algorithm  consists  of  two  phases.  The  first  (Phase  I)  uses  an  exter- 
nal squared  penalty  function  to  find  a point  close  to  a local  minimum. 

u 


Starting  with  the  algorithm  then  solves  a sequence  of  linearly  constrain- 

u 

ed  problems  (Phase  II).  Selected  nonlinear  constraints  are  linearized  for 
each  such  Phase  II  Iteration.  With  suitable  assumptions,  convergence  from 
any  initial  point,  with  quadratic  convergence  in  Phase  II,  is  shown. 

The  practical  implementation  of  this  algorithm  is  described,  and  its 
potential  application  to  a model  for  the  assessment  of  energy  alternatives 
is  discussed  briefly. 


UNCLASSIFIED ^ 

■tCUNITV  CLAtUPICATlOM  THK  ^AairVkM  OM« 


